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We discuss applications of the theory of Quantum Chaos to one of the paradigm models 
of many-body quantum physics — the Bose-Hubbard model, which describes, in par¬ 
ticular, interacting ultracold Bose atoms in an optical lattice. After preliminary, pure 
quantum analysis of the system we introduce the classical counterpart of the Bose- 
Hubbard model and the governing semiclassical equations of motion. We analyze these 
equations for the problem of Bloch oscillations of cold atoms where a number of exper¬ 
imental results are available. The review is written for non-experts and can be viewed 
as an introduction to the field. 
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I. INTRODUCTION 


The Bose-Hubbard (BH) model, introduced by Gersch and Knollman in 1963 (Gersch and Knollman 1963), became 
a hot topic in physics after the theoretical work by Jaksch et. al. of 1998 (Jaksch et al. |1998 ), where it was noticed 
that this model can be realized with the help of cold atoms loaded into an optical lattice, and, especially, after the 


fundamental experiment by Greiner et. al. of 2002 (Greiner et ail 2002), where the authors demonstrated quantum 


phase transition from the super-fluid state to the Mott insulator for rubidium atoms in a cubic optical lattice. Ever 


since the BH model is discussed almost exclusively in the context of cold atoms physics (Jaksch and Zoller 2005). 


In the past decade many different phenomena of the cold atoms physics, which are described by the Bose-Hubbard 
and Bose-Hubbard like Hamiltonians, were studied in great detail. Naturally, these studies contributed to our under¬ 
standing of the properties of the model, which is currently considered as paradigm model of many-body physics that 
can be tested experimentally. Still, we are far from complete understanding and every novel analytical approach or 
laboratory experiment adds something new to our knowledge of the BH model. 

In this tutorial review we describe a new approach to the BH model which is based on the ideas of Quantum Ghaos 
- a modern theory which deals with quantum non-integrable systems (Giannoni et al. 1991 Stockmann, 1999). A 


particular feature of this approach is the extensive use of classical mechanics. This might come as a surprise because 
the BH model is usually considered as a genuine quantum system, with no classical counterpart. However, as it will 
be explained in the review, there is an analogy between the quantum and classical descriptions of a single-particle 
system, and the microscopic and mean-field descriptions of a many-body system. This analogy helps us to better 
understand the BH model, especially, when it concerns excited states of the system. 

The review consists of three parts. In Sec. [H] we consider the BH model as a generic complex quantum system, 
without appealing to the classical mechanics. Following the main idea of Quantum Ghaos that the energy spectrum 
of a complex system should have common features with the spectrum of random matrices, we perform statistical 
analysis of eigenvalues and eigenfunctions of the BH Hamiltonian and compare the result with RMT (Random Matrix 
Theory) predictions. Section HI begins with discussing of the semiclassical limit, where we follow the method of the 
truncated Husimi function. The power of this semiclassical method is demonstrated in Sec. |IV.A[ where we derive 
the Bogoliubov spectrum by using semiclassical quantization, and in Sec. |IV[ where we consider the problem of Bloch 
oscillations (BOs) of cold atoms in tilted ID optical lattices. It is shown, in particular, that one can reproduce 
quantum dynamics of the system by solving classical equations. 
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II. ENERGY SPECTRUM OF THE BOSE-HUBBARD MODEL 
A. The Bose-Hubbard Hamiltonian 

Having in mind cold Bose atoms in a ID optical lattice the BH Hamiltonian reads 


Hn = -i 


(4+iaj + /i-c.) + ^ ^ n; (n; - 1) . 

Z=1 Z=1 


( 1 ) 


In Eq. Q the index I labels the lattice sites (wells of the optical potential), o; and dj are the bosonic annihilation 
and creation operators, 


fii = djdi is the number operator, J the hopping matrix element (the rate of inter-well tunneling), and U the 
microscopic interaction constant (the energy paid by two atoms sharing the same well). The constant U is mainly 
determined by the s-wave scattering length as for neutral atoms, 


U = 


Airagd? 

M 


and the constant J by the lattice depth Vq, 


I </'/(!■) I r 




J = J (l>i+i{z)Hs(j)i{z)dz , Hs = - +VoCOs'^{kLz) 

(here (t>i{z) are the Wannier functions localized at Hh well of the optical potential and fci, is the laser wave vector). 
Both the scattering length and the depth Vq can be varied in large intervals, which allows practically arbitrary ratio 
U/J. 



z/d 

FIG. 1 Pictorial presentation of cold atoms (open circles) in an optical lattice. The thin green line shows one of the Wannier 
function cj)i{z). 


^ One-dimensional optical lattice or, more exactly, an array of independent ID lattices is created by using two strong standing laser waves 
in the x and y directions, which create the so-called quantum tubes, and one weak standing wave in the z direction, which periodically 
modulates the quantum tubes. 









4 





•\ 






\ 

\ 


\ ^ 

\. *'v 




VU-.., 






\ 

\ 

\ 





\ \ 





\\> 



\ 




\%\ 


\ 



\ 

\ 



\ ^ 



\ ^ 






FIG. 2 Example of the Hamiltonian matrix for N = L = 5. Shown are nonzero matrix elements of the Hamiltonians 0 and 
(j^ in the Fock basis Q and (|^, respectively. 


A remark concerning boundary conditions is in turn. In a laboratory experiment the default boundary condition is 
residual harmonic confinement due to finite widths of the laser beams. In the theory, however, one usually considers 
the periodic boundary condition, where (L + l)th site of the lattice is identified with the first site, i.e., al+i = 0 

Through the paper we shall assume the later case. Notice that the periodic boundary condition implies conservation 
of the total quasimomentum. This can be seen by rewriting the Hamiltonian Q in terms of the operators bk and , 


h 1 \" f 2TTk 


ai 


h' - 

"fc — 




( 2 ) 


Operators 0 annihilate or create an atom in the Bloch state with the quasimomentum k = 2'iTk/L. Using the 
transformation 0 we have 

Ho = -JY^ cos ^ XI + ^2 - fca - ki) , (3) 

k '' ' ki,k2,k3,k4 


where S is the periodic (^-function, i.e., d{k) equals unity if A: is a multiple of L and zero otherwise. The presence of 
the 5-function in the interaction term insures that the total quasimomentum is conserved. 

Finally, we discuss the Hilbert space of the BH system. It is spanned by the Fock states 

|n) = |ni,n 2 ,... ,nL) , = N , (4) 

i 


where N is the total number of atoms. The dimension of the Hilbert space is 


M 


[N + L-l)\ 
N\{L-1)\ 


In the coordinate representation the basis state @ is given by the symmetrized product of N Wannier functions 
4>i{z). Correspondently, if we consider the basis state of the Hamiltonian 0, 

|n) = |ni,n 2 ,... ,71^) , '^Uk = N , (5) 

k 


2 


We mention, in passing, that few-site BH model with periodic boundary condition can be realized experimentally by using non-trivial 
Gaussian beams llAmico et aL||2005|l. This setup, however, excludes consideration of the limit L ^ oo. 
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FIG. 3 Energy levels of the BH system {N = L = 5) as functions of the parameter u where J = 1 — u and U = u . The left 
panel shows the whole spect rum, the right panel is independ ent subset of odd states with zero total quasimomentum. The 
figure is borrowed from Ref. (Kolovsky and Buchleitn'^ 20041. 


it is given by the symmetrized product of N Bloch waves ^k{z) = L ^/‘^'^iexp{i2TTkl/L)(j)i{z). The total quasimo¬ 
mentum of the state © is calculated as 

K = ^modi (6) 

and can take one of L values. Knowing the action of bosonic operators on a given Fock state, 

= y/ni\ ..., n; - 1,...) , ajl ... ,ni,...) = vW+Tl ... ,ni + 1,...) , 

(for operators bk and bl we have similar equations) we calculate the Hamiltonian matrix of the size Af x Af. The 
explicit form of this matrix depends on particular ordering of the basis states. The main point, however, is that this 
matrix appears to be very sparse, see Fig. ([^. 

B. Statistical analysis of the energy spectrum 

First we give a numerical evidence that the BH model belongs to the class of quantum non-integrable systems. 
Figure shows the energy spectrum of the system for L — N = 5 which is parametrized by the parameter u, where 
J = 1 — u and U = u. This spectrum can be decomposed into L independent spectra associated with different values 
of the total quasimomentum. In the Bloch representation ra one finds independent spectra by considering different 
subsets of the Hilbert space labeled by the quasimomentum . Alternatively, we can find these spectra by calculating 
the matrix of the Hamiltonian Q in the translationally invariant basis 

L 

|m, k) = {1/v^)^ exp(m/)5'^|ni,n2,. -. ,ni) , 

1=1 

where k is the total quasimomentum and S is the cyclic permutation operator: S\ni,n 2 , ■ ■ ■ ,ni) = |n 2 ,n 3 ,... 

The spectrum associated with k = 0 can be decomposed further according to the reflection (odd-even) symmetry. 
At this stage the decomposition is complete. As an example, the right panel in Fig. shows one of the independent 
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FIG. 4 Upper panel: The mean density of states, histogram, approximated by Eq. 0- solid line. Lower panel: Integrated 
distribution I{s) = P{s')ds' for one of independent spectra as compared to the Poisson and Wigner-Dyson distributions. 
The system parameters are L = N = 8 and u = 0.3. 


spectra. It is seen that the energy levels exhibit avoided crossings as they approach each other. This proves that the 
BH model is not integrable. Moreover, it appears to be a chaotic system in the sense of Quantum Chaos. To prove 
the last statement one should perform statistical analysis of the spectrum, where the simplest test is the distribution 
of normalized distances between the nearest levels. 

To normalize the distances or, using the Quantum Chaos terminology, to unfold the spectrum we need to know the 
mean density of states p{E). An estimate for p{E) can be obtained as follows. Let us set for the moment U = 0. 
Then the spectrum is known analytically [see Eq. ([^]: 



Here are integer numbers and, if > 1, this should be viewed as the sum of equal terms. Thus the energy 
level E is given by a sum of N real numbers, each restricted by modulus by J. Considering the thermodynamic limit 
N,L ^ oo and using the central limit theorem we conclude that levels E are distributed according to the normal law 
with the variance ~ JN. The normal distribution proves to be a good approximation also for nonzero U < J. In 
this case the Gaussian is shifted as the whole by the mean interaction energy: 


p{E) - exp 


(E — Eint) 


21 


^JN , Er, 


UN^ 


(7) 


Validity of the approximation Q is illustrated in the upper panel in Fig. A more thorough analysis reveals 
deviations of the actual density of states from Eq. Q, especially at the tails of the distribution. However, for the sake 
of statistical analysis of the spectrum, where the main contribution comes from the central part of the spectrum, it 
is quite satisfactory. 

Having the mean density of state obtained we introduce the normalized distance between the nearest energy levels. 


s — {En-\-i Eji') pi^Eyi) , 


( 8 ) 


® Exclusions are the cases C/ = 0 or J = 0 and the case of two-site BH model (L = 2), where the spectrum can be found analytically for 
arbitrary U and J. 
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FIG. 5 Panels (a) and (b): Density of states for N = 7, L ='3, J=l, e = 0.1, and the interaction constant 
U = 0.2, respectivel y. Panel ( c): In tegrated level spacing distributions for the central part of the spectrum, 
borrowed from Ref. ( Kolovsk^ 20061. 


U = 0.02 and 
The figure is 


calculate the distribution function P{s) for the distance Q, and compare it with the Wigner-Dyson distributions for 
random matrices. This comparison reveals remarkable agreement with P(s) for the Gaussian Orthogonal Ensemble 
(GOE) of random matrices, 

P(s) = |sexp , (9) 

see lower panel in Fig. Thus the system 0 is chaotic in the sense of Quantum Chaos. It will be shown later on in 
Sec. |III| that the BH model is a chaotic system also in the sense of Classical Chaos. This explains the amazing fact 
that the sparse matrix shown in Fig. has similar properties as a fully random matrix. 


C. Transition to chaos 


The transition to chaos in the BH model takes place as a transition over the parameter U/J, which alone defines 
the properties of the system. To study this transition we can avoid complex procedure of the spectrum decomposition 
by introducing a weak on-site disorder, 


H = Ho 


L 

E 

z=i 


eini 


|d| < e 


( 10 ) 


If the disorder is weak enough, so that the Anderson localization length is much larger than the system size L, the 
main effect of the disorder is destruction of the global symmetries (which are the total quasimomentum and odd-even 
symmetry). Thus we may study the spectrum statistics without preliminary decomposition of the spectrum. 

Figure illustrates transition to chaos in the system (10) as we increase the interaction constant U. One clearly 
observes the change from the Poisson statistics. 


P{s) = exp(-s) , 


( 11 ) 


These are real symmetric matrices with random entries according to the normal low. More exactly, probability density to meet matrix 
H in the ensemble is given by V{H) ~ ex.p[—ATr{H^)] where A is the normalization constant. A common choice is A = 7r^/2A/* where 
JV is the matrix size. Then the mean density of states of the GOE matrix is given by p{E) = y/l — {7rE/2Af)'^. 
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FIG. 6 Gray-scale image of the matrix (12|. The system parameters are the same as in Fig. 


which is typical for integrable systems, to the Wigner-Dyson statistics ([^, which is a hallmark of quantum chaotic 


systems (Stockmann 19991. 


Another indication of the transition to chaos comes from analysis of the eigenstates Note that the only 

quantum number of the eigenstate |'I'„) is its energy and we assume that the states are ordered according to their 
energies, i.e., En > Em if n > m. A useful characteristic of eigenstates is the matrix 

i?(m,n) = |(T'^(17')|^n(17))|2 , (12) 

where U and U' are two different values of our control parameter. Clearly, R{n, m) is the identity matrix ii U = U'. 
However, if U' deviates from U it become a banded matrix, see Fig. The crucial point is that for chaotic systems 
i?(n, m) obeys the universal distribution 

r/27r 

(n — m)2 -h r^/d ’ 

1996). In this equation the parameter F is a function of the 


R{n — m) = 


(13) 


known as the Breit-Wigner equation (Fyodorov et al. 


difference U' — U and the bar denotes an average over several eigenstates. (Without this averaging procedure, the 
matrix elements R{n,m) show strong fluctuations, see Fig. [^) 

The Breit-Wigner distribution (13) is illustrated in Fig. In this figure the dots are numerical data and the solid 
line is Eq. (13). A reasonable agreement is noticed. We mention that numerical data in Fig. are collected from the 
central part of the spectrum which is marked by the solid lines in Fig.[^a,b). For the low- and high-energy eigenstates 
one finds essential deviations from Eq. (13). This observation tells that these states are not chaotic or, at least, not 


fully chaotic. We come back to this point in Sec. |III.C| 

To conclude this section we briefly discuss a transition from GOE to GUE (Gaussian Unitary Ensemble) [^spectrum 
statistics in the BH model. The GOE statistics is typical for quantum systems with time-reversal symmetry, which is 
obviously the case of the Hamiltonian 0- One can break this symmetry by introducing the complex hopping matrix 
element: 


H = -i 
2 


Lj 


'i+i 


aie 


,tS 


i=i 


U 


h-c.) E-^fiiifii - 1) . 


(14) 




^ This matrix is closely related to the so-called local density of states, R{m, E) = i?(m, n)5[E — En), which has a number of important 

physical applications, see Ref. ||Kolovsky| |2009|| , for example. 

® This is the ensemble of hermitian randorn matrices with the probability density 'P(H) ~ ex.p[—ATr{H'^H)]. 
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FIG. 7 Mean values of the matrix elements across the main diagonal in the central part of the matrix on the linear (left panel) 
and logarithmic (right panel) scales. The solid line is the best fit by the Breit-Wigner formula. 


The Hamiltonian (14) appears in the problem of atomic Bloch oscillations and can be experimentally realized by 


shaking the lattice with a proper frequency ( 

Kolovsky 

r 

>011 

I. The effect of non-zero phase 9 on dynamics of cold 

atoms was studied, for example, in Ref. ( 

Haller et al. 

201( 

11. In this section our prime interest is the energy spectrum, 


p. ^ 32 2 / 4 2 

P{s) = -^s exp- s 

TT^ \ TT 


as 9 deviates from zero. This, however, does not affect the main result of this section that the BH model is a quantum 
chaotic system. 


III. SEMICLASSICAL QUANTIZATION OF THE BOSE-HUBBARD MODEL 

In the theory of Quantum Chaos analysis of a quantum system is usually preceded by the analysis of its classical 
counterpart. In this review we reverted this sequence for the reason that the quantum analysis is actually simpler. 
Now we come to the classical consideration, where the first step is to identify the classical counterpart of the BH 
model. 


A. Semiclassical limit 

There are several ways to introduce classical counterpart of the quantum system 0 - We shall follow the approach 
based on the notion of the Husimi function. Below we shall introduce the effective Planck constant heff which is 
inverse proportional to number of atoms, hgff = 1/N. This constant should not be mismatched with the fundamental 
Planck constant h which we set to unity from now on. 


^ A similar approach is based on the notion of the Wigner function ( |Polkovnikov[ |2010[ [Sinatra et al.\ |2002| [Steel et al.\ [l998[ |. The 
Husimi function, however, has an advantage that it is positively defined. 


























10 


Given |'I'(i)) to be the many-body wave function of the quantum Hamiltonian, the Husimi function is defined as 

/(a,t) = |(a|^h(f))^ (15) 

where |a) are the so-called coherent SU{L) states ( Perelomo^ 1986), 

L 


1 




E' 


N 


aia\ I \vac) . 


(15) the Schrodinger equation for the wave function |4'(t)) takes the form 

df / 1 


Note that the Husimi function (|15|) is a function of L complex amplitude a/ and time. In terms of the Husimi function 

(16) 

(17) 


where denotes the Poisson brackets, the c-number Hamiltonian Hq reads 


Ho = -^ ^( 4 + 10 ; -f c.c.) + f XI’ 9= ^ , 


1=1 


and we refer the reader to the work ( 

Trimborn et al. 

2008 

) for explicit form of the terms which are inverse proportional 

to N. The constant g in the classical Hamiltonian ( 

17 

) is called the macroscopic interaction constant, to distinguish 


it from the microscopic interaction constant U. 

We remark that, formally, the Hamiltonian 0 follows from the microscopic Hamiltonian ([^ by using the ‘quan¬ 
tization rules’ ai -f-)- a* -O- and H o i7/h, where n = N/L is the filling factor (the mean number 

of atoms per lattice site). The approach of the Husimi function (as well as the approach of the Wigner function) 
provides a justification of these quantization rules. 

Let us now consider the mean-field limit N —)■ cx), C/ —)■ 0 while g = const. In this limit the terms proportional to 
1/N vanish and the Husimi function reduces to the multidimensional ^-function. 


/(a,t) = 


1=1 


where ai{t) satisfies the Hamilton equations of motion. 


d dtdiQ J , \ I I I 

= vnr = (fli-i-i + ai-i) + g\ai\ 


dt 


daf 


'ai 


(18) 


Eq. (181 is known in the physical literature as the Discrete Nonlinear Schrodinger Equation (DNLSE) and can be 


viewed as a discrete analogue of the Gross-Pitaevskii equation for a Bose-Einstein condensate (Pitaevskii and Stringari 


2003). We mention that DNLSE also appears in problems of Nonlinear Optics and molecular vibrations, where it has 


been studied for decades (Eilbeck et al. 1985). 


A comment is due on the terminology. In the title of this section we used ‘semiclassical limit’ instead of the ‘mean- 


field limit’. The reason is that Eq. (16) formally coincides with equation on the Husimi function of a single-particle 


system if one identifies 1/A^ with the Planck constant. Thus one can use the common semiclassical theory to study 
the BH model, - we shall give an example in Sec. |HLG| 

Finally, we mention that within the formalism of the Husimi function we are actually not bound with the limit 
N ^ oo and may consider finite N as well. Of course, Eq. (16) is simpler than the original Schrodinger equation for 


the many-body wave function |'I'(t)) only if we neglect the last term in this equation - the approximation known as the 
truncated Husimi function. Fortunately, there are physically important cases where this approximation is justified. 
In these cases the truncated Husimi function correctly reproduces quantum dynamics of the BH model even when the 


mean-field equation (18) fails to do this. 


B. Phase space of the classical Bose-Hubbard system 


Eq. (18) defines classical trajectories in 2L dimensional phase space which lie on the energy shell E = i7(a). 
Depending on the energy E and initial conditions the trajectory can be either regular or chaotic. (Here we don’t 
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discuss the trivial case L = 2 where all trajectories are regular.) It is largely an open question to hnd the volume 
of regular and chaotic components on a given energy shell. In the general case of L-site system we only know that 
there are regular islands for E close to the energy of the ground. As it will be shown in the next subsection, these 
islands are associated with the Bogoliubov spectrum for elementary excitations of a Bose-Einstein condensate (BEC). 
For higher energy shells, even if there are stability islands, their size is expected to be small. Statistical analysis of 
the energy spectrum of the quantum BH model presented in Sec. [Tl] strongly supports this statement. In fact, if there 
were large regular islands for higher energies, we would see this as a deviation in the spectrum statistics from the 
Wigner-Dyson distribution. 


C. Bogoliubov spectrum 


Let us discuss the low-energy stability islands in more detail. To make analysis simpler we shall consider L = 3. 
Generalization to larger L is straighforward and consists of substituting the quantity 6 = J[I — cos(27r/3)] in the 
equations below by the quantity 6k = T[I — cos(27rfc/L)]. 

The analysis involves several steps. First we rewrite the Hamiltonian 0 for L = 3 in terms of the canonical 
variables bk and b^. This gives 


Hr 


= -^E 


cos 




/ 27rA:\ 


Kh + ^ 


- y 

2 ^ 

ki,k2,ks,k4 


-I- A:2 - ^3 - ki) . 


Next we switch to the action-angle variables, bk = y/7k e^pii4>k), and explicitly take into account that J2k^k = 1- 
This reduces our system of three degrees of freedom to a system of two degrees of freedom: 

Hq = (<5 -f g)(I-i + I+i) + 2 ^/ 0 "s/d-i/+i cos(^_i + 4>+i) (19) 

+ I^i + /+i) + 2g cos{2(f>zfi — (j)±i) , 

± 

where 6 = J[I — cos(27r/3)], Iq = 1 — /_i — I+i and the phases (j)±i of variables b±i{t) are measured with respect to 
the phase of bo{t). The low-energy dynamics of the system (191, which is associated with the low-energy spectrum of 


the quantum system, implies I±i <C Iq- Keeping in the Hamiltonian (19) only the terms linear in I±i, and using one 
more canonical transformation, 

I = I+i + I-i , 0 = {(jj+i + y’-i)/2 , 

M = (J_|_i — /_i)/2 , 0 = (p+i — 4>-i , 

we obtain the effective Hamiltonian which locally describes the low-energy stability island: 

H^ff = {5 + g)I + g\/p - 4M2 cos(26») , \M\<I/2. (20) 

Note that H^ff does not include phase 0 and, hence, the action M is a constant of motion. 

The obtained Hamiltonian (20) suffices to find the low-energy spectrum of the 3-site BH model. To do this we 


integrate the system (20) by introducing new action, I = (l/27r) f I{9, E)d9, and resolving this equation with respect 
to the energy. We get 

E = ni , n = yj2g6 + (52 , (21) 

where the frequency H is nothing else as the Bogoliubov frequency. Finally, we quantize actions I and M in units 
of the effective Planck constant fieff = l/A^- This gives equidistant set of energy levels En = Eq + iln, with (n -|- 1) 
degeneracy of the nth level. 

It is interesting to compare the above result with the exact energy spectrum of the 3-site model. This spectrum is 
shown in Fig. where, to facilitate the comparison, we rescale it by using the Bogoliubov frequency fl. As expected, 
the degenerate equidistant spectrum is a good approximation only up to some critical energy, above which the classical 
dynamics of the BH system is chaotic. In principle, one can use the semiclassical quantization also in the chaotic 
region. However, this will require more sophisticated semiclassical theory known as the Periodic Orbits Theory [see 


Chapter 7 in Ref. (Stockmann 1999)], which is based on the notion of the van Vleck-Gutzwiller propagator. An 
application of the van Vleck-Gutzwiller propagator to the BH model is found in Ref. (Engl et al. 2014). 


® More results are known for the 3-site system, where the classical phase space can be relatively easy visualized, see recent work ( |Arwas| 
I et aZ.j |2014[ | and references therein. 

® For L > 3 the system has several Bogoliubov frequencies Qk = \J‘igSk + <5^ ~ y/gk, where k is usually interpreted as the wave vector of 
elementary excitations. 
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FIG. 8 Energy spectrum of the 3-site BH model for N — 40. The energy is measured with respect to the ground energy Eq 
and scaled with respect to the Bogoliubov frequency fl. Error bars indicate energy intervals, where the classical counterpart of 
the system shows chaotic dynamic. The figure is borrowed from Ref. ([Kolovsl^ |2007[). 


IV. BLOCH OSCILLATIONS OF BOSE ATOMS 


In this section we discuss Bloch oscillations (BOs) of interacting Bose atoms. The main reason for discussing this 
phenomenon in the present review is that BOs can test chaotic nature of the BH model. Furthermore, recently BOs 
of Bose atoms in ID lattices have been studied experimentally (Meinert et al. 20141, providing the first experimental 
results on Quantum Chaos in the BH model. 


A. Governing equation 

In the single-band approximation (which is assumed throughout the paper) BOs of interacting atoms are described 
by the Hamiltonian 

H = HoEdF^lni, ( 22 ) 

I 

where Hq is the BH Hamiltonian Q, F a static (for example, gravitational) field, and d the lattice period (d = 1 in 
what follows). Introducing the Bloch frequency ujb = dF/Ti = F and using the substitution di —> diexp{—iFlt) the 
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time-independent Hamiltonian (22) reduces to the form (14) where 9 = Ft: 

L 


~ 2 ^ (“i-i-i 

1=1 


aLi h.c^ +^^ni{hi - 1) . 

^ 1=1 


(23) 


Referring to the experiment (Meinert et al. [2014 ) the initial wave function of interacting atoms is given by the ground 
state of the Hamiltonian Hq = H{t = 0) and the simplest quantity to be measured is the mean momentum per one 
atom: 

J 


pit) = ■ 


Notice that for non-interacting atoms we would have p{t) = — J sin(i^t). This equation is the essence of phenomenon 
of BOs. The problem to be addressed is the effect of atom-atom interactions which, as it was shown in Sec. [TJ make 
the Hamiltonian (23) chaotic. 


B. Mean-field analysis 

Let us first analyze BOs of interacting atoms by using the mean-field approach: 


= -^ (ai+ie*^‘ -b a/_ie -b g\ai\'^ai , 


It is easy to check that Eq. (24) has a periodic solution E3 


ai(t) = exp ( i— sm{Ft) — igt 


(24) 


(25) 


For the mean momentum per atom Eq. (25) gives 


pit) = -j}^ 


= — Jsin(Et) . 


analysis of the periodic trajectory (25). 


Thus one can consider the solution (25) as a candidate for BOs of interacting atoms, where the next step is stability 


The stability analysis is done in the usual way, i.e., by linearizing Eq. (24) around the periodic trajectory (25): 


z —5a = Al[a(t)]5a 


(26) 


(here 6a{t) is a deviation from the periodic trajectory and A4 the Jacobi matrix). It leads to the following result 


(Kolovsky et al. 2009 Zheng et al. 2004). In the limit of large L the parameter space of the system (24) is divided 


into two parts by the critical line 


p 

CA 


( 3g , F <2J 
{ y/TOgJ , F>2J ■ 


(27) 


In the strong field regime F > F^r all Lyapunov exponents of the linear Eq. (26) are zero and, hence, the solution 


(25) is stable. In the opposite case F < F^r there are positive exponents and the solution is unstable. Pictorially, 


this result is illustrated in Fig. which shows periodic trajectory in the multi-dimensional phase-space of the system 
together with a nearby trajectory. 


Not periodic phase exp(—igt) in Eq. l|25ll is irrelevant and can be removed by the obvious substitution. 
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FIG. 9 Pictorial presentation of the stable and unstable periodic trajectory in the multi-dimensional phase-space of the classical 
BH system. 


C. Quantum ensemble 


We have shown that the trajectory (25), which is defined by the initial conditions ai{t = 0) = 1, is unstable in 
the weak field regime. Thus an arbitrary small deviation from the specified initial conditions results in a completely 
different trajectory. This brings us back to the Husimi function, which describes the evolution of not a single 
trajectory but of an ensemble of trajectories. This ensemble is obviously defined by the equation 


/(a,t = 0) = |(a|4'(t = 0))|^ 


(28) 


where |'I'(t = 0)) is the initial many-body wave function given by the ground state of the Hamiltonian Hq. For U < J 
and h ~ 1 this ground state is well approximated by a BEG of non-interacting atoms, 

= 0)) = \vac) , ^ a} . (29) 


Then the distribution (281 is known analytically and we can generate ensemble of initial conditions by using, for 


example, the exception-rejection numerical method (Press et al.. 2007). In what follows we shall refer to this ensemble 


of initial conditions as the quantum ensemble, to stress that it is defined by the quantum many-body state of the 


system. As an example. Fig. 10 shows the quantum ensemble for the state (29) in the Wannier representations for 


L = 5 and N = 15. It is seen that amplitudes ai{t = 0) deviates from unity by both the absolute value and the phase, 
where the characteristic size of deviations is inverse proportional to -s/h- 


Having the quantum ensemble in hands we evolve each trajectory according to Eq. (24) and average the result over 


the ensemble. (Clearly, this amounts to the Monte-Carlo solution of the truncated equation on the Husimi function.) 
The panel (a) in Fig. [Tl| shows the mean atomic momentum as function of time for the static field F < Fcr- An 
exponential decay of BOs, 


p{t) = — J exp(—yt) sin(F't) 


(30) 


is clearly seen. The panel (a) should be compared with the panel (b) showing the solution of the Schrodinger equation 
with the microscopic Hamiltonian (23). An excellent agreement is noticed. We discuss the physics behind this 


remarkable result in the next subsection. 


D. Internal decoherence 


Previous studies of the Bloch dynamics of interacting Bose atoms (]Buchleitner and Kolovsky 2003 Kolovsky et 


al. 2009) proved that exponential decay of BOs is due to decoherence of the initial BEC state or self-thermalization 


11 


This phenomenon is often referred to as the dynamical or modulation instability. 
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phase/it 

FIG. 10 Quantum ensemble (100 realizations) representing the many-body BEG state for = 15 and L = 5. 
widths of distributions are inverse proportional to 


The characteristic 





FIG. 11 Bloch oscillations of interacting atoms calculated by using the truncated Husimi function (left) and the full quantum- 
mechanical simulations (right). The system parameters are N = 15, L = 5, J = 1, U = 0.1/3, and F = 0.1 (top), and E = 10 
(bottom). 


of the system. On the formal level the self-thermalization means that the one-particle density matrix TZ{t) relaxes to 
a diagonal (in the Bloch basis) matrix with equal populations of the single-particle quasimomentum modes. Thus the 
linear entropy of the system S = Tr(7^^), which is one of possible characteristics of the system coherence, decreases 
from unity to S' = 1/L ^ 1. From the viewpoint of classical mechanics the self-thermalization is a consequence of 
chaotic dynamics of the system or, more precisely, the mixing property of the chaotic dynamics. 

The discussed ‘internal’ decoherence has common features with ‘external’ decoherence caused by interaction of the 
system with the environment. In particular, if we consider BOs of a single atom coupled to a bath, 0 we shall also 


12 


In the problem of atomic BOs the bath is given by zero modes of the electromagnetic field which are responsible for spontaneous 
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observe the exponential decay of oscillations (Kolovsky et al. 2002). On the other hand, we know from studies on the 


general problem of quantum-classical correspondence that external decoherence suppresses interference terms in the 
governing equation of motion (Kolovsky 1994 1996 Zurek 1991|). In Eq. (16) on the Husimi function these terms 


are denoted as 0{1/N). The excellent agreement between exact quantum simulations and the semiclassical approach 
of the truncated Husimi function proves that this is also the case for the internal decoherence. Thus we have a loop: 
classical chaotic dynamics is responsible for the internal decoherence which, in its turn, causes the quantum system 
to behave classically. 

It is interesting to consider the strong field limit F > Fcr which breaks the above loop. Now the classical dynamics 
is stable and, hence, there is no internal decoherence. Dynamics of the mean momentum for F > F^r is shown in 
the panel (d) in Fig. 11 It presents periodic revivals of BOs which are described by the following simple equation 
dKolovskyl [2^ , 


p{t) = —J exp (—2n[l — cos{Ut)]) sin{Ft) , n = N/L . 


(31) 


It is seen in Fig. 113") that the truncated Husimi function approach does not reproduce the revivals, which are a 


quantum interference effect. Recent laboratory studies of BOs in ID lattices (Meinert et al. 2014) undoubtedly 


confirm the transition from the quasi-periodic dynamics (31) to the exponentially decaying oscillations (30) as the 
static field F is decreased. F3 


V. CONCLUSIONS 


We discussed the Bose-Hubbard Hamiltonian beyond its ground state. It was proven that this system is a chaotic 
system in the sense of Quantum Chaos. Of course, it is not the only many-body system which is chaotic - the other 


examples are provided, for instance, by non-integrable models of spin chains (Berman et al. 2001 Montambaux et 


al. 1993). However, the BH model is extremely important in the experimental cold-atom physics and for this reason 


deserves a special attention. 

In the review we focused on the case where the kinetic energy of atoms dominates the interaction energy {J > U). In 
this case the ground state of the system is known to be a super-fluid state with low-energy excitations described by the 
Bogoliubov theory. We revisited this problem from uncommon perspective of the quantum-classical correspondence 
and showed that the Bogoliubov spectrum can be obtained by quantizing the low-energy stability islands of the 
classical BH model. 

If we go to higher energies, the energy spectrum of the BH model become irregular and must be analyzed statistically. 
We considered the simplest statistical characteristic of the spectrum - the level spacing distribution - which was shown 
to obey the Wigner-Dyson equation for random matrices. 

Another direction of research is response of the BH system to a sudden change of its parameters (so-called quench 
dynamics) or to external perturbations. In the review we considered the response to a static field. For cold atoms in 
optical lattices this could be, for example, the gravitational field (Poll et al.. 2011). For the initial condition given by 


the ground state of the system, the static held induces Bloch oscillations which, according to the mean-held analysis, 
can be either stable or unstable. The microscopic analysis of the system shows that the unstable regime results in the 
exponential decay of BOs due to internal decoherence, which is a consequence of the chaotic dynamics of the system. 
Remarkably, in this case we were able to reproduce the quantum dynamics by using the ‘classical’ approach of the 
truncated Husimi function. This result sheds the new light on the old problem of the onset of classicality in our world. 

We conclude this review by mentioning some future prospects. As stated above, we focussed on the case J > U 
where the ground state of the BH model is a super-fluid state. It is interesting to study the opposite case U ^ J (where 
the ground state of the system is a Mott insulator) by using a semiclassical approach. In particular, this concerns 
the problem of quantum phase transition from the Mott-insulator state to the density-wave state in tilted ID lattices 
(Sachdev et al.. 2002 Simon et al. 2011). The second avenue is a generalizing the results to two-dimensional case 


where, besides potential fields, one can include into consideration gauge fields (Aidelsburger et al. 2011,2013 Miyake 


et al. 2013). Since atoms are neutral these fields are often referred to as synthetic magnetic fields. For vanishing 


interactions dynamical and spectral properties of cold atoms in a 2D lattice subject to the synthetic magnetic and 


electric fields are discussed in detail in the recent work (Kolovsky and Mantica 2014). 


emission. In laboratory experiments intensity of this process and, hence, the rate of external decoherence is controlled by tuning the 
laser frequency further or closer to the atomic resonance. 

In the cited experiment the authors used gravitational field which was compensated by the levitation force to a desired level. 
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